﻿using System;
using System.Collections.Generic;
using System.Linq;
using System.Text;
using System.Threading.Tasks;
using SpiceNet.Circuits;
using SpiceNet.Diagnostics;
using System.Numerics;
using SpiceNet.Parameters;

namespace SpiceNet.Models
{
    public class BSIM1 : Instance
    {
        /// <summary>
        /// Get parameters
        /// </summary>
        public override Dictionary<string, IP> ParamTable => _pTable;

        #region Parameters
        private enum _c
        {
            /* device parameters */
            BSIM1_W = 1,
            BSIM1_L = 2,
            BSIM1_AS = 3,
            BSIM1_AD = 4,
            BSIM1_PS = 5,
            BSIM1_PD = 6,
            BSIM1_NRS = 7,
            BSIM1_NRD = 8,
            BSIM1_OFF = 9,
            BSIM1_IC_VBS = 10,
            BSIM1_IC_VDS = 11,
            BSIM1_IC_VGS = 12,
            BSIM1_IC = 13,
            /* device questions */
            BSIM1_DNODE = 201,
            BSIM1_GNODE = 202,
            BSIM1_SNODE = 203,
            BSIM1_BNODE = 204,
            BSIM1_DNODEPRIME = 205,
            BSIM1_SNODEPRIME = 206,
            BSIM1_VBD = 207,
            BSIM1_VBS = 208,
            BSIM1_VGS = 209,
            BSIM1_VDS = 210,
            BSIM1_CD = 211,
            BSIM1_CBS = 212,
            BSIM1_CBD = 213,
            BSIM1_GM = 214,
            BSIM1_GDS = 215,
            BSIM1_GMBS = 216,
            BSIM1_GBD = 217,
            BSIM1_GBS = 218,
            BSIM1_QB = 219,
            BSIM1_CQB = 220,
            BSIM1_QG = 221,
            BSIM1_CQG = 222,
            BSIM1_QD = 223,
            BSIM1_CQD = 224,
            BSIM1_CGG = 225,
            BSIM1_CGD = 226,
            BSIM1_CGS = 227,
            BSIM1_CBG = 228,
            BSIM1_CAPBD = 231,
            BSIM1_CQBD = 232,
            BSIM1_CAPBS = 233,
            BSIM1_CQBS = 234,
            BSIM1_CDG = 235,
            BSIM1_CDD = 236,
            BSIM1_CDS = 237,
            BSIM1_VON = 238,
            BSIM1_QBS = 239,
            BSIM1_QBD = 240,
            BSIM1_SOURCECONDUCT = 241,
            BSIM1_DRAINCONDUCT = 242
        }
        private static Dictionary<string, IP> _pTable = new Dictionary<string, IP>
        {
            { "l", new IP(IF.IOP, (int)_c.BSIM1_L, T.REAL, "Length") },
            { "w", new IP(IF.IOP, (int)_c.BSIM1_W, T.REAL, "Width") },
            { "ad", new IP(IF.IOP, (int)_c.BSIM1_AD, T.REAL, "Drain area") },
            { "as", new IP(IF.IOP, (int)_c.BSIM1_AS, T.REAL, "Source area") },
            { "pd", new IP(IF.IOP, (int)_c.BSIM1_PD, T.REAL, "Drain perimeter") },
            { "ps", new IP(IF.IOP, (int)_c.BSIM1_PS, T.REAL, "Source perimeter") },
            { "nrd", new IP(IF.IOP, (int)_c.BSIM1_NRD, T.REAL, "Number of squares in drain") },
            { "nrs", new IP(IF.IOP, (int)_c.BSIM1_NRS, T.REAL, "Number of squares in source") },
            { "off", new IP(IF.IOP, (int)_c.BSIM1_OFF, T.FLAG, "Device is initially off") },
            { "vds", new IP(IF.IOP, (int)_c.BSIM1_IC_VDS, T.REAL, "Initial D-S voltage") },
            { "vgs", new IP(IF.IOP, (int)_c.BSIM1_IC_VGS, T.REAL, "Initial G-S voltage") },
            { "vbs", new IP(IF.IOP, (int)_c.BSIM1_IC_VBS, T.REAL, "Initial B-S voltage") },
            { "ic", new IP(IF.IP, (int)_c.BSIM1_IC, T.VECTOR, "Vector of DS,GS,BS initial voltages") }
        };

        public Parameter<double> B1vdsat { get; } = new Parameter<double>();
        public Parameter<double> B1von { get; } = new Parameter<double>();

        public Parameter<double> B1l { get; } = new Parameter<double>(5e-6);
        public Parameter<double> B1w { get; } = new Parameter<double>(5e-6);
        public Parameter<double> B1drainArea { get; } = new Parameter<double>();
        public Parameter<double> B1sourceArea { get; } = new Parameter<double>();
        public Parameter<double> B1drainPerimeter { get; } = new Parameter<double>();
        public Parameter<double> B1sourcePerimeter { get; } = new Parameter<double>();
        public Parameter<double> B1drainSquares { get; } = new Parameter<double>(1.0);
        public Parameter<double> B1sourceSquares { get; } = new Parameter<double>(1.0);
        public bool B1off { get; set; } = false;
        public Parameter<double> B1icVDS { get; } = new Parameter<double>();
        public Parameter<double> B1icVGS { get; } = new Parameter<double>();
        public Parameter<double> B1icVBS { get; } = new Parameter<double>();

        // Normal variables
        public double B1GDoverlapCap, B1GSoverlapCap, B1GBoverlapCap, B1drainConductance, B1sourceConductance, B1vfb, B1phi, B1K1, B1K2, B1eta, B1etaB, B1etaD, B1betaZero, B1betaZeroB, B1ugs, B1ugsB, B1uds,
            B1udsB, B1udsD, B1betaVdd, B1betaVddB, B1betaVddD, B1subthSlope, B1subthSlopeB, B1subthSlopeD, B1vt0;
        #endregion

        /// <summary>
        /// Instance equation indices
        /// </summary>
        private int B1dNode;
        private int B1gNode;
        private int B1sNode;
        private int B1bNode;
        private int B1dNodePrime;
        private int B1sNodePrime;
        private int B1state;
        private int B1mode = 1;

        #region Constants
        private const int B1vbd = 0;
        private const int B1vbs = 1;
        private const int B1vgs = 2;
        private const int B1vds = 3;
        private const int B1cd = 4;
        private const int B1id = 4;
        private const int B1cbs = 5;
        private const int B1ibs = 5;
        private const int B1cbd = 6;
        private const int B1ibd = 6;
        private const int B1gm = 7;
        private const int B1gds = 8;
        private const int B1gmbs = 9;
        private const int B1gbd = 10;
        private const int B1gbs = 11;
        private const int B1qb = 12;
        private const int B1cqb = 13;
        private const int B1iqb = 13;
        private const int B1qg = 14;
        private const int B1cqg = 15;
        private const int B1iqg = 15;
        private const int B1qd = 16;
        private const int B1cqd = 17;
        private const int B1iqd = 17;
        private const int B1cggb = 18;
        private const int B1cgdb = 19;
        private const int B1cgsb = 20;
        private const int B1cbgb = 21;
        private const int B1cbdb = 22;
        private const int B1cbsb = 23;
        private const int B1capbd = 24;
        private const int B1iqbd = 25;
        private const int B1cqbd = 25;
        private const int B1capbs = 26;
        private const int B1iqbs = 27;
        private const int B1cqbs = 27;
        private const int B1cdgb = 28;
        private const int B1cddb = 29;
        private const int B1cdsb = 30;
        private const int B1vono = 31;
        private const int B1vdsato = 32;
        private const int B1qbs = 33;
        private const int B1qbd = 34;
        #endregion

        /// <summary>
        /// Constructor
        /// </summary>
        /// <param name="name"></param>
        public BSIM1(string name) : base(name, 4, typeof(BSIM1Model)) { }

        /// <summary>
        /// Set parameter
        /// </summary>
        protected override void Param(int id, object value, Circuit ckt = null)
        {
            double[] v;
            switch ((_c)id)
            {
                case _c.BSIM1_W:
                    B1w.Par((double)value);
                    break;
                case _c.BSIM1_L:
                    B1l.Par((double)value);
                    break;
                case _c.BSIM1_AS:
                    B1sourceArea.Par((double)value);
                    break;
                case _c.BSIM1_AD:
                    B1drainArea.Par((double)value);
                    break;
                case _c.BSIM1_PS:
                    B1sourcePerimeter.Par((double)value);
                    break;
                case _c.BSIM1_PD:
                    B1drainPerimeter.Par((double)value);
                    break;
                case _c.BSIM1_NRS:
                    B1sourceSquares.Par((double)value);
                    break;
                case _c.BSIM1_NRD:
                    B1drainSquares.Par((double)value);
                    break;
                case _c.BSIM1_OFF:
                    B1off = (bool)value;
                    break;
                case _c.BSIM1_IC_VBS:
                    B1icVBS.Par((double)value);
                    break;
                case _c.BSIM1_IC_VDS:
                    B1icVDS.Par((double)value);
                    break;
                case _c.BSIM1_IC_VGS:
                    B1icVGS.Par((double)value);
                    break;
                case _c.BSIM1_IC:
                    v = (double[])value;
                    switch (v.Length)
                    {
                        case 3:
                            B1icVBS.Par(v[2]); goto case 2;
                        case 2:
                            B1icVGS.Par(v[1]); goto case 1;
                        case 1:
                            B1icVDS.Par(v[0]); break;
                        default:
                            throw new BadParameterException();
                    }
                    break;
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Ask parameter
        /// </summary>
        protected override object Ask(int id, Circuit ckt = null)
        {
            switch ((_c)id)
            {
                case _c.BSIM1_L:
                    return B1l;
                case _c.BSIM1_W:
                    return B1w;
                case _c.BSIM1_AS:
                    return B1sourceArea;
                case _c.BSIM1_AD:
                    return B1drainArea;
                case _c.BSIM1_PS:
                    return B1sourcePerimeter;
                case _c.BSIM1_PD:
                    return B1drainPerimeter;
                case _c.BSIM1_NRS:
                    return B1sourceSquares;
                case _c.BSIM1_NRD:
                    return B1drainSquares;
                case _c.BSIM1_OFF:
                    return B1off;
                case _c.BSIM1_IC_VBS:
                    return B1icVBS;
                case _c.BSIM1_IC_VDS:
                    return B1icVDS;
                case _c.BSIM1_IC_VGS:
                    return B1icVGS;
                case _c.BSIM1_DNODE:
                    return B1dNode;
                case _c.BSIM1_GNODE:
                    return B1gNode;
                case _c.BSIM1_SNODE:
                    return B1sNode;
                case _c.BSIM1_BNODE:
                    return B1bNode;
                case _c.BSIM1_DNODEPRIME:
                    return B1dNodePrime;
                case _c.BSIM1_SNODEPRIME:
                    return B1sNodePrime;
                case _c.BSIM1_SOURCECONDUCT:
                    return B1sourceConductance;
                case _c.BSIM1_DRAINCONDUCT:
                    return B1drainConductance;
                case _c.BSIM1_VBD:
                    return ckt.State.States[0][B1state + B1vbd];
                case _c.BSIM1_VBS:
                    return ckt.State.States[0][B1state + B1vbs];
                case _c.BSIM1_VGS:
                    return ckt.State.States[0][B1state + B1vgs];
                case _c.BSIM1_VDS:
                    return ckt.State.States[0][B1state + B1vds];
                case _c.BSIM1_CD:
                    return ckt.State.States[0][B1state + B1cd];
                case _c.BSIM1_CBS:
                    return ckt.State.States[0][B1state + B1cbs];
                case _c.BSIM1_CBD:
                    return ckt.State.States[0][B1state + B1cbd];
                case _c.BSIM1_GM:
                    return ckt.State.States[0][B1state + B1gm];
                case _c.BSIM1_GDS:
                    return ckt.State.States[0][B1state + B1gds];
                case _c.BSIM1_GMBS:
                    return ckt.State.States[0][B1state + B1gmbs];
                case _c.BSIM1_GBD:
                    return ckt.State.States[0][B1state + B1gbd];
                case _c.BSIM1_GBS:
                    return ckt.State.States[0][B1state + B1gbs];
                case _c.BSIM1_QB:
                    return ckt.State.States[0][B1state + B1qb];
                case _c.BSIM1_CQB:
                    return ckt.State.States[0][B1state + B1cqb];
                case _c.BSIM1_QG:
                    return ckt.State.States[0][B1state + B1qg];
                case _c.BSIM1_CQG:
                    return ckt.State.States[0][B1state + B1cqg];
                case _c.BSIM1_QD:
                    return ckt.State.States[0][B1state + B1qd];
                case _c.BSIM1_CQD:
                    return ckt.State.States[0][B1state + B1cqd];
                case _c.BSIM1_CGG:
                    return ckt.State.States[0][B1state + B1cggb];
                case _c.BSIM1_CGD:
                    return ckt.State.States[0][B1state + B1cgdb];
                case _c.BSIM1_CGS:
                    return ckt.State.States[0][B1state + B1cgsb];
                case _c.BSIM1_CBG:
                    return ckt.State.States[0][B1state + B1cbgb];
                case _c.BSIM1_CAPBD:
                    return ckt.State.States[0][B1state + B1capbd];
                case _c.BSIM1_CQBD:
                    return ckt.State.States[0][B1state + B1cqbd];
                case _c.BSIM1_CAPBS:
                    return ckt.State.States[0][B1state + B1capbs];
                case _c.BSIM1_CQBS:
                    return ckt.State.States[0][B1state + B1cqbs];
                case _c.BSIM1_CDG:
                    return ckt.State.States[0][B1state + B1cdgb];
                case _c.BSIM1_CDD:
                    return ckt.State.States[0][B1state + B1cddb];
                case _c.BSIM1_CDS:
                    return ckt.State.States[0][B1state + B1cdsb];
                case _c.BSIM1_VON:
                    return ckt.State.States[0][B1state + B1vono];
                case _c.BSIM1_QBS:
                    return ckt.State.States[0][B1state + B1qbs];
                case _c.BSIM1_QBD:
                    return ckt.State.States[0][B1state + B1qbd];
                default:
                    throw new BadParameterException();
            }
        }

        /// <summary>
        /// Setup the device
        /// </summary>
        /// <param name="ckt"></param>
        public override void Setup(Circuit ckt, Model mod, object p = null)
        {
            BSIM1Model model = mod as BSIM1Model;

            // Create extra nodes if the contacts have a conductance
            List<NodeType> extra = new List<NodeType>();

            bool addDrain = false;
            if ((model.B1sheetResistance != 0.0) && (this.B1drainSquares != 0.0))
            {
                extra.Add(NodeType.Voltage);
                addDrain = true;
            }

            bool addSource = false;
            if ((model.B1sheetResistance != 0.0) && (this.B1drainSquares != 0.0))
            {
                extra.Add(NodeType.Voltage);
                addSource = true;
            }

            // Map nodes
            var nodes = BindNodes(ckt, extra.ToArray());
            B1dNode = nodes[0].Number;
            B1gNode = nodes[1].Number;
            B1sNode = nodes[2].Number;
            B1bNode = nodes[3].Number;
            int index = 4;
            B1dNodePrime = addDrain ? nodes[index++].Number : B1dNode;
            B1sNodePrime = addSource ? nodes[index++].Number : B1sNode;

            B1state = CreateState(ckt, 35);

        }

        /// <summary>
        /// Do temperature dependent
        /// </summary>
        /// <param name="ckt"></param>
        public override void Temperature(Circuit ckt, Model mod, object p = null)
        {
            BSIM1Model model = mod as BSIM1Model;
            double EffChanLength;
            double EffChanWidth;
            double CoxWoverL;
            double Leff;    /* effective channel length im micron */
            double Weff;    /* effective channel width in micron */


            if ((EffChanLength = this.B1l - model.B1deltaL * 1e-6) <= 0)
                throw new CircuitException(String.Format("B1: mosfet {0}, model {1}: Effective channel length <=0", Name, model.Name));
            if ((EffChanWidth = this.B1w - model.B1deltaW * 1e-6) <= 0)
                throw new CircuitException(String.Format("B1: mosfet {0}, model {1}: Effective channel width <=0", Name, model.Name));
            this.B1GDoverlapCap = EffChanWidth * model.B1gateDrainOverlapCap;
            this.B1GSoverlapCap = EffChanWidth * model.B1gateSourceOverlapCap;
            this.B1GBoverlapCap = this.B1l * model.B1gateBulkOverlapCap;

            /* process drain series resistance */
            if ((this.B1drainConductance = model.B1sheetResistance *
                    this.B1drainSquares) != 0.0)
            {
                this.B1drainConductance = 1.0 / this.B1drainConductance;
            }

            /* process source series resistance */
            if ((this.B1sourceConductance = model.B1sheetResistance *
                    this.B1sourceSquares) != 0.0)
            {
                this.B1sourceConductance = 1.0 / this.B1sourceConductance;
            }

            Leff = EffChanLength * 1.0e6; /* convert into micron */
            Weff = EffChanWidth * 1.0e6; /* convert into micron */
            CoxWoverL = model.B1Cox * Weff / Leff; /* F/cm**2 */

            this.B1vfb = model.B1vfb0 +
                model.B1vfbL / Leff + model.B1vfbW / Weff;
            this.B1phi = model.B1phi0 +
                model.B1phiL / Leff + model.B1phiW / Weff;
            this.B1K1 = model.B1K10 +
                model.B1K1L / Leff + model.B1K1W / Weff;
            this.B1K2 = model.B1K20 +
                model.B1K2L / Leff + model.B1K2W / Weff;
            this.B1eta = model.B1eta0 +
                model.B1etaL / Leff + model.B1etaW / Weff;
            this.B1etaB = model.B1etaB0 +
                model.B1etaBl / Leff + model.B1etaBw / Weff;
            this.B1etaD = model.B1etaD0 +
                model.B1etaDl / Leff + model.B1etaDw / Weff;
            this.B1betaZero = model.B1mobZero;
            this.B1betaZeroB = model.B1mobZeroB0 +
                model.B1mobZeroBl / Leff + model.B1mobZeroBw / Weff;
            this.B1ugs = model.B1ugs0 +
                model.B1ugsL / Leff + model.B1ugsW / Weff;
            this.B1ugsB = model.B1ugsB0 +
                model.B1ugsBL / Leff + model.B1ugsBW / Weff;
            this.B1uds = model.B1uds0 +
                model.B1udsL / Leff + model.B1udsW / Weff;
            this.B1udsB = model.B1udsB0 +
                model.B1udsBL / Leff + model.B1udsBW / Weff;
            this.B1udsD = model.B1udsD0 +
                model.B1udsDL / Leff + model.B1udsDW / Weff;
            this.B1betaVdd = model.B1mobVdd0 +
                model.B1mobVddl / Leff + model.B1mobVddw / Weff;
            this.B1betaVddB = model.B1mobVddB0 +
                model.B1mobVddBl / Leff + model.B1mobVddBw / Weff;
            this.B1betaVddD = model.B1mobVddD0 +
                model.B1mobVddDl / Leff + model.B1mobVddDw / Weff;
            this.B1subthSlope = model.B1subthSlope0 +
                model.B1subthSlopeL / Leff + model.B1subthSlopeW / Weff;
            this.B1subthSlopeB = model.B1subthSlopeB0 +
                model.B1subthSlopeBL / Leff + model.B1subthSlopeBW / Weff;
            this.B1subthSlopeD = model.B1subthSlopeD0 +
                model.B1subthSlopeDL / Leff + model.B1subthSlopeDW / Weff;

            if (this.B1phi < 0.1) this.B1phi = 0.1;
            if (this.B1K1 < 0.0) this.B1K1 = 0.0;
            if (this.B1K2 < 0.0) this.B1K2 = 0.0;

            this.B1vt0 = this.B1vfb + this.B1phi + this.B1K1 *
                Math.Sqrt(this.B1phi) - this.B1K2 * this.B1phi;

            this.B1von.Value = this.B1vt0;  /* added for initialization*/ // NOTE: Overwrites model parameter...?

            /* process Beta Parameters (unit: A/V**2) */

            this.B1betaZero = this.B1betaZero * CoxWoverL;
            this.B1betaZeroB = this.B1betaZeroB * CoxWoverL;
            this.B1betaVdd = this.B1betaVdd * CoxWoverL;
            this.B1betaVddB = this.B1betaVddB * CoxWoverL;
            this.B1betaVddD = Math.Max(this.B1betaVddD * CoxWoverL, 0.0);
        }

        /// <summary>
        /// Load the BSIM1 model
        /// </summary>
        /// <param name="ckt"></param>
        public override void Load(Circuit ckt, Model mod, object p = null)
        {
            BSIM1Model model = mod as BSIM1Model;
            CircuitState state = ckt.State;
            Integration.IntegrationMethod integ = ckt.Integration;

            double DrainSatCurrent = 0.0;
            double EffectiveLength = 0.0;
            double GateBulkOverlapCap = 0.0;
            double GateDrainOverlapCap = 0.0;
            double GateSourceOverlapCap = 0.0;
            double SourceSatCurrent = 0.0;
            double DrainArea = 0.0;
            double SourceArea = 0.0;
            double DrainPerimeter = 0.0;
            double SourcePerimeter = 0.0;
            double arg = 0.0;
            double capbd = 0.0;
            double capbs = 0.0;
            double cbd = 0.0;
            double cbhat = 0.0;
            double cbs = 0.0;
            double cd = 0.0;
            double cdrain = 0.0;
            double cdhat = 0.0;
            double cdreq = 0.0;
            // double ceq = 0.0;
            double ceqbd = 0.0;
            double ceqbs = 0.0;
            double ceqqb = 0.0;
            double ceqqd = 0.0;
            double ceqqg = 0.0;
            double czbd = 0.0;
            double czbdsw = 0.0;
            double czbs = 0.0;
            double czbssw = 0.0;
            double delvbd = 0.0;
            double delvbs = 0.0;
            double delvds = 0.0;
            double delvgd = 0.0;
            double delvgs = 0.0;
            double evbd = 0.0;
            double evbs = 0.0;
            double gbd = 0.0;
            double gbs = 0.0;
            double gcbdb = 0.0;
            double gcbgb = 0.0;
            double gcbsb = 0.0;
            double gcddb = 0.0;
            double gcdgb = 0.0;
            double gcdsb = 0.0;
            double gcgdb = 0.0;
            double gcggb = 0.0;
            double gcgsb = 0.0;
            double gcsdb = 0.0;
            double gcsgb = 0.0;
            double gcssb = 0.0;
            double gds = 0.0;
            // double geq = 0.0;
            double gm = 0.0;
            double gmbs = 0.0;
            double sarg = 0.0;
            double sargsw = 0.0;
            // double tol = 0.0;
            double vbd = 0.0;
            double vbs = 0.0;
            double vcrit = 0.0;
            double vds = 0.0;
            double vdsat = 0.0;
            double vgb = 0.0;
            double vgd = 0.0;
            double vgdo = 0.0;
            double vgs = 0.0;
            double von = 0.0;
            // double xfact = 0.0;
            double xnrm = 0.0;
            double xrev = 0.0;
            bool Check = false;
            double cgdb = 0.0;
            double cgsb = 0.0;
            double cbdb = 0.0;
            double cdgb = 0.0;
            double cddb = 0.0;
            double cdsb = 0.0;
            double cggb = 0.0;
            double cbgb = 0.0;
            double cbsb = 0.0;
            double csgb = 0.0;
            double cssb = 0.0;
            double csdb = 0.0;
            double PhiB = 0.0;
            double PhiBSW = 0.0;
            double MJ = 0.0;
            double MJSW = 0.0;
            double argsw = 0.0;
            double qgate = 0.0;
            double qbulk = 0.0;
            double qdrn = 0.0;
            double qsrc = 0.0;
            double cqgate = 0.0;
            double cqbulk = 0.0;
            double cqdrn = 0.0;
            double vt0 = 0.0;
            double[] args = new double[8];
            int ByPass;

            EffectiveLength = B1l - model.B1deltaL * 1.0e-6;/* m */
            DrainArea = B1drainArea;
            SourceArea = B1sourceArea;
            DrainPerimeter = B1drainPerimeter;
            SourcePerimeter = B1sourcePerimeter;
            if ((DrainSatCurrent = DrainArea * model.B1jctSatCurDensity)
                    < 1e-15)
            {
                DrainSatCurrent = 1.0e-15;
            }
            if ((SourceSatCurrent = SourceArea * model.B1jctSatCurDensity)
                    < 1.0e-15)
            {
                SourceSatCurrent = 1.0e-15;
            }
            GateSourceOverlapCap = model.B1gateSourceOverlapCap * B1w;
            GateDrainOverlapCap = model.B1gateDrainOverlapCap * B1w;
            GateBulkOverlapCap = model.B1gateBulkOverlapCap * EffectiveLength;
            von = model.B1type * B1von;
            vdsat = model.B1type * B1vdsat;
            vt0 = model.B1type * B1vt0;

            Check = true;
            ByPass = 0;
            if (((ckt.Mode & Circuit.Modes.InitSmSig) != 0))
            {
                vbs = state.States[0][B1state + B1vbs];
                vgs = state.States[0][B1state + B1vgs];
                vds = state.States[0][B1state + B1vds];
            }
            else if (((ckt.Mode & Circuit.Modes.InitTran) != 0))
            {
                vbs = state.States[1][B1state + B1vbs];
                vgs = state.States[1][B1state + B1vgs];
                vds = state.States[1][B1state + B1vds];
            }
            else if (((ckt.Mode & Circuit.Modes.InitJct) != 0) && !B1off)
            {
                vds = model.B1type * B1icVDS;
                vgs = model.B1type * B1icVGS;
                vbs = model.B1type * B1icVBS;
                if ((vds == 0) && (vgs == 0) && (vbs == 0) &&
                        (((ckt.Mode &
                        (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.DcOp | Circuit.Modes.DcTranCurve)) != 0) ||
                        (!((ckt.Mode & Circuit.Modes.Uic) != 0))))
                {
                    vbs = -1;
                    vgs = vt0;
                    vds = 0;
                }
            }
            else if (((ckt.Mode & (Circuit.Modes.InitJct | Circuit.Modes.InitFix)) != 0) &&
                  (B1off))
            {
                vbs = vgs = vds = 0;
            }
            else
            {
                vbs = model.B1type * (
                    state.Solution[B1bNode] -
                    state.Solution[B1sNodePrime]);
                vgs = model.B1type * (
                    state.Solution[B1gNode] -
                    state.Solution[B1sNodePrime]);
                vds = model.B1type * (
                    state.Solution[B1dNodePrime] -
                    state.Solution[B1sNodePrime]);
                vbd = vbs - vds;
                vgd = vgs - vds;
                vgdo = state.States[0][B1state + B1vgs] -
                    state.States[0][B1state + B1vds];
                delvbs = vbs - state.States[0][B1state + B1vbs];
                delvbd = vbd - state.States[0][B1state + B1vbd];
                delvgs = vgs - state.States[0][B1state + B1vgs];
                delvds = vds - state.States[0][B1state + B1vds];
                delvgd = vgd - vgdo;

                if (B1mode >= 0)
                {
                    cdhat =
                        state.States[0][B1state + B1cd] -
                        state.States[0][B1state + B1gbd] * delvbd +
                        state.States[0][B1state + B1gmbs] * delvbs +
                        state.States[0][B1state + B1gm] * delvgs +
                        state.States[0][B1state + B1gds] * delvds;
                }
                else
                {
                    cdhat =
                        state.States[0][B1state + B1cd] -
                        (state.States[0][B1state + B1gbd] -
                          state.States[0][B1state + B1gmbs]) * delvbd -
                        state.States[0][B1state + B1gm] * delvgd +
                        state.States[0][B1state + B1gds] * delvds;
                }
                cbhat =
                    state.States[0][B1state + B1cbs] +
                    state.States[0][B1state + B1cbd] +
                    state.States[0][B1state + B1gbd] * delvbd +
                    state.States[0][B1state + B1gbs] * delvbs;

                von = model.B1type * B1von;
                if (state.States[0][B1state + B1vds] >= 0)
                {
                    vgs = fetlim(vgs, state.States[0][B1state + B1vgs]
                            , von);
                    vds = vgs - vgd;
                    vds = limvds(vds, state.States[0][B1state + B1vds]);
                    vgd = vgs - vds;
                }
                else
                {
                    vgd = fetlim(vgd, vgdo, von);
                    vds = vgs - vgd;
                    vds = -limvds(-vds, -(state.States[0][B1state + B1vds]));
                    vgs = vgd + vds;
                }
                if (vds >= 0)
                {
                    vcrit = Circuit.CONSTvt0 * Math.Log(Circuit.CONSTvt0 / (Circuit.CONSTroot2 * SourceSatCurrent));
                    vbs = pnjlim(vbs, state.States[0][B1state + B1vbs],
                            Circuit.CONSTvt0, vcrit, ref Check); /* B1 test */
                    vbd = vbs - vds;
                }
                else
                {
                    vcrit = Circuit.CONSTvt0 * Math.Log(Circuit.CONSTvt0 / (Circuit.CONSTroot2 * DrainSatCurrent));
                    vbd = pnjlim(vbd, state.States[0][B1state + B1vbd],
                            Circuit.CONSTvt0, vcrit, ref Check); /* B1 test*/
                    vbs = vbd + vds;
                }
            }

            /* determine DC current and derivatives */
            vbd = vbs - vds;
            vgd = vgs - vds;
            vgb = vgs - vbs;


            if (vbs <= 0.0)
            {
                gbs = SourceSatCurrent / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbs = gbs * vbs;
            }
            else
            {
                evbs = Math.Exp(vbs / Circuit.CONSTvt0);
                gbs = SourceSatCurrent * evbs / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbs = SourceSatCurrent * (evbs - 1) + ckt.Config.Gmin * vbs;
            }
            if (vbd <= 0.0)
            {
                gbd = DrainSatCurrent / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbd = gbd * vbd;
            }
            else
            {
                evbd = Math.Exp(vbd / Circuit.CONSTvt0);
                gbd = DrainSatCurrent * evbd / Circuit.CONSTvt0 + ckt.Config.Gmin;
                cbd = DrainSatCurrent * (evbd - 1) + ckt.Config.Gmin * vbd;
            }
            /* line 400 */
            if (vds >= 0)
            {
                /* normal mode */
                B1mode = 1;
            }
            else
            {
                /* inverse mode */
                B1mode = -1;
            }
            /* call B1evaluate to calculate drain current and its 
             * derivatives and charge and capacitances related to gate
             * drain, and bulk
             */
            if (vds >= 0)
            {
                B1evaluate(ckt, model, vds, vbs, vgs, ref gm, ref gds, ref gmbs, ref qgate,
                    ref qbulk, ref qdrn, ref cggb, ref cgdb, ref cgsb, ref cbgb, ref cbdb, ref cbsb, ref cdgb,
                    ref cddb, ref cdsb, ref cdrain, ref von, ref vdsat);
            }
            else
            {
                B1evaluate(ckt, model, -vds, vbd, vgd, ref gm, ref gds, ref gmbs, ref qgate,
                    ref qbulk, ref qsrc, ref cggb, ref cgsb, ref cgdb, ref cbgb, ref cbsb, ref cbdb, ref csgb,
                    ref cssb, ref csdb, ref cdrain, ref von, ref vdsat);
            }

            B1von.Value = model.B1type * von;
            B1vdsat.Value = model.B1type * vdsat;



            /*
             *  COMPUTE EQUIVALENT DRAIN CURRENT SOURCE
             */
            cd = B1mode * cdrain - cbd;
            if (((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac | Circuit.Modes.InitSmSig)) != 0) ||
                    (((ckt.Mode & Circuit.Modes.TranOp) != 0) &&
                    ((ckt.Mode & Circuit.Modes.Uic) != 0)))
            {
                /*
                 *  charge storage elements
                 *
                 *   bulk-drain and bulk-source depletion capacitances
                 *  czbd : zero bias drain junction capacitance
                 *  czbs : zero bias source junction capacitance
                 * czbdsw:zero bias drain junction sidewall capacitance
                 * czbssw:zero bias source junction sidewall capacitance
                 */

                czbd = model.B1unitAreaJctCap * DrainArea;
                czbs = model.B1unitAreaJctCap * SourceArea;
                czbdsw = model.B1unitLengthSidewallJctCap * DrainPerimeter;
                czbssw = model.B1unitLengthSidewallJctCap * SourcePerimeter;
                PhiB = model.B1bulkJctPotential;
                PhiBSW = model.B1sidewallJctPotential;
                MJ = model.B1bulkJctBotGradingCoeff;
                MJSW = model.B1bulkJctSideGradingCoeff;

                /* Source Bulk Junction */
                if (vbs < 0)
                {
                    arg = 1 - vbs / PhiB;
                    argsw = 1 - vbs / PhiBSW;
                    sarg = Math.Exp(-MJ * Math.Log(arg));
                    sargsw = Math.Exp(-MJSW * Math.Log(argsw));
                    state.States[0][B1state + B1qbs] =
                        PhiB * czbs * (1 - arg * sarg) / (1 - MJ) + PhiBSW *
                    czbssw * (1 - argsw * sargsw) / (1 - MJSW);
                    capbs = czbs * sarg + czbssw * sargsw;
                }
                else
                {
                    state.States[0][B1state + B1qbs] =
                        vbs * (czbs + czbssw) + vbs * vbs * (czbs * MJ * 0.5 / PhiB
                        + czbssw * MJSW * 0.5 / PhiBSW);
                    capbs = czbs + czbssw + vbs * (czbs * MJ / PhiB +
                        czbssw * MJSW / PhiBSW);
                }

                /* Drain Bulk Junction */
                if (vbd < 0)
                {
                    arg = 1 - vbd / PhiB;
                    argsw = 1 - vbd / PhiBSW;
                    sarg = Math.Exp(-MJ * Math.Log(arg));
                    sargsw = Math.Exp(-MJSW * Math.Log(argsw));
                    state.States[0][B1state + B1qbd] =
                        PhiB * czbd * (1 - arg * sarg) / (1 - MJ) + PhiBSW *
                    czbdsw * (1 - argsw * sargsw) / (1 - MJSW);
                    capbd = czbd * sarg + czbdsw * sargsw;
                }
                else
                {
                    state.States[0][B1state + B1qbd] =
                        vbd * (czbd + czbdsw) + vbd * vbd * (czbd * MJ * 0.5 / PhiB
                        + czbdsw * MJSW * 0.5 / PhiBSW);
                    capbd = czbd + czbdsw + vbd * (czbd * MJ / PhiB +
                        czbdsw * MJSW / PhiBSW);
                }

            }

            /*
             *  check convergence
             */
            if ((!B1off) || (!((ckt.Mode & Circuit.Modes.InitFix) != 0)))
            {
                if (Check)
                {
                    ckt.IsCon = false;
                }
            }
            state.States[0][B1state + B1vbs] = vbs;
            state.States[0][B1state + B1vbd] = vbd;
            state.States[0][B1state + B1vgs] = vgs;
            state.States[0][B1state + B1vds] = vds;
            state.States[0][B1state + B1cd] = cd;
            state.States[0][B1state + B1cbs] = cbs;
            state.States[0][B1state + B1cbd] = cbd;
            state.States[0][B1state + B1gm] = gm;
            state.States[0][B1state + B1gds] = gds;
            state.States[0][B1state + B1gmbs] = gmbs;
            state.States[0][B1state + B1gbd] = gbd;
            state.States[0][B1state + B1gbs] = gbs;

            state.States[0][B1state + B1cggb] = cggb;
            state.States[0][B1state + B1cgdb] = cgdb;
            state.States[0][B1state + B1cgsb] = cgsb;

            state.States[0][B1state + B1cbgb] = cbgb;
            state.States[0][B1state + B1cbdb] = cbdb;
            state.States[0][B1state + B1cbsb] = cbsb;

            state.States[0][B1state + B1cdgb] = cdgb;
            state.States[0][B1state + B1cddb] = cddb;
            state.States[0][B1state + B1cdsb] = cdsb;

            state.States[0][B1state + B1capbs] = capbs;
            state.States[0][B1state + B1capbd] = capbd;

            /* bulk and channel charge plus overlaps */

            if ((!((ckt.Mode & (Circuit.Modes.Tran | Circuit.Modes.Ac)) != 0)) &&
                    ((!((ckt.Mode & Circuit.Modes.TranOp) != 0)) ||
                    (!((ckt.Mode & Circuit.Modes.Uic) != 0))) && (!((ckt.Mode
                    & Circuit.Modes.InitSmSig) != 0))) goto line850;

            // line755:
            if (B1mode > 0)
            {

                args[0] = GateDrainOverlapCap;
                args[1] = GateSourceOverlapCap;
                args[2] = GateBulkOverlapCap;
                args[3] = capbd;
                args[4] = capbs;
                args[5] = cggb;
                args[6] = cgdb;
                args[7] = cgsb;

                B1mosCap(ckt, vgd, vgs, vgb,
            args,
                    /*
                    GateDrainOverlapCap,
                    GateSourceOverlapCap,GateBulkOverlapCap,
                    capbd,capbs,
                                cggb,cgdb,cgsb,
                    */
            cbgb, cbdb, cbsb,
            cdgb, cddb, cdsb,
                        ref gcggb, ref gcgdb, ref gcgsb,
            ref gcbgb, ref gcbdb, ref gcbsb,
            ref gcdgb, ref gcddb, ref gcdsb, ref gcsgb, ref gcsdb, ref gcssb,
            ref qgate, ref qbulk,
                        ref qdrn, ref qsrc);
            }
            else
            {

                args[0] = GateSourceOverlapCap;
                args[1] = GateDrainOverlapCap;
                args[2] = GateBulkOverlapCap;
                args[3] = capbs;
                args[4] = capbd;
                args[5] = cggb;
                args[6] = cgsb;
                args[7] = cgdb;

                B1mosCap(ckt, vgs, vgd, vgb,
            args,
                    /*
                    GateSourceOverlapCap,
                    GateDrainOverlapCap,GateBulkOverlapCap,
                    capbs,capbd,
                    cggb,cgsb,cgdb,
                    */
            cbgb, cbsb, cbdb,
            csgb, cssb, csdb,
            ref gcggb, ref gcgsb, ref gcgdb,
            ref gcbgb, ref gcbsb, ref gcbdb,
            ref gcsgb, ref gcssb, ref gcsdb, ref gcdgb, ref gcdsb, ref gcddb,
            ref qgate, ref qbulk,
            ref qsrc, ref qdrn);
            }

            if (ByPass != 0) goto line860;
            state.States[0][B1state + B1qg] = qgate;
            state.States[0][B1state + B1qd] = qdrn -
                    state.States[0][B1state + B1qbd];
            state.States[0][B1state + B1qb] = qbulk +
                    state.States[0][B1state + B1qbd] +
                    state.States[0][B1state + B1qbs];

            /* store small signal parameters */
            if ((!((ckt.Mode & (Circuit.Modes.Ac | Circuit.Modes.Tran)) != 0)) &&
                    ((ckt.Mode & Circuit.Modes.TranOp) != 0) && ((ckt.Mode &
                    Circuit.Modes.Uic) != 0)) goto line850;
            if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
            {
                state.States[0][B1state + B1cggb] = cggb;
                state.States[0][B1state + B1cgdb] = cgdb;
                state.States[0][B1state + B1cgsb] = cgsb;
                state.States[0][B1state + B1cbgb] = cbgb;
                state.States[0][B1state + B1cbdb] = cbdb;
                state.States[0][B1state + B1cbsb] = cbsb;
                state.States[0][B1state + B1cdgb] = cdgb;
                state.States[0][B1state + B1cddb] = cddb;
                state.States[0][B1state + B1cdsb] = cdsb;
                state.States[0][B1state + B1capbd] = capbd;
                state.States[0][B1state + B1capbs] = capbs;

                goto line1000;
            }

            if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                state.States[1][B1state + B1qb] =
                    state.States[0][B1state + B1qb];
                state.States[1][B1state + B1qg] =
                    state.States[0][B1state + B1qg];
                state.States[1][B1state + B1qd] =
                    state.States[0][B1state + B1qd];
            }

            integ.Integrate(state, B1state + B1qb, 0.0);
            integ.Integrate(state, B1state + B1qg, 0.0);
            integ.Integrate(state, B1state + B1qd, 0.0);

            goto line860;

        line850:
            /* initialize to zero charge conductance and current */
            ceqqg = ceqqb = ceqqd = 0.0;
            gcdgb = gcddb = gcdsb = 0.0;
            gcsgb = gcsdb = gcssb = 0.0;
            gcggb = gcgdb = gcgsb = 0.0;
            gcbgb = gcbdb = gcbsb = 0.0;
            goto line900;

        line860:
            /* evaluate equivalent charge current */
            cqgate = state.States[0][B1state + B1iqg];
            cqbulk = state.States[0][B1state + B1iqb];
            cqdrn = state.States[0][B1state + B1iqd];
            ceqqg = cqgate - gcggb * vgb + gcgdb * vbd + gcgsb * vbs;
            ceqqb = cqbulk - gcbgb * vgb + gcbdb * vbd + gcbsb * vbs;
            ceqqd = cqdrn - gcdgb * vgb + gcddb * vbd + gcdsb * vbs;

            if ((ckt.Mode & Circuit.Modes.InitTran) != 0)
            {
                state.States[1][B1state + B1iqb] =
                    state.States[0][B1state + B1iqb];
                state.States[1][B1state + B1iqg] =
                    state.States[0][B1state + B1iqg];
                state.States[1][B1state + B1iqd] =
                    state.States[0][B1state + B1iqd];
            }

            /*
             *  load current vector
             */
        line900:

            ceqbs = model.B1type * (cbs - (gbs - ckt.Config.Gmin) * vbs);
            ceqbd = model.B1type * (cbd - (gbd - ckt.Config.Gmin) * vbd);

            ceqqg = model.B1type * ceqqg;
            ceqqb = model.B1type * ceqqb;
            ceqqd = model.B1type * ceqqd;
            if (B1mode >= 0)
            {
                xnrm = 1;
                xrev = 0;
                cdreq = model.B1type * (cdrain - gds * vds - gm * vgs - gmbs * vbs);
            }
            else
            {
                xnrm = 0;
                xrev = 1;
                cdreq = -(model.B1type) * (cdrain + gds * vds - gm * vgd - gmbs * vbd);
            }

            state.Rhs[B1gNode] -= ceqqg;
            state.Rhs[B1bNode] -= (ceqbs + ceqbd + ceqqb);
            state.Rhs[B1dNodePrime] +=
                    (ceqbd - cdreq - ceqqd);
            state.Rhs[B1sNodePrime] +=
                    (cdreq + ceqbs + ceqqg + ceqqb + ceqqd);

            /*
             *  load y matrix
             */

            state.Matrix[B1dNode, B1dNode] += (B1drainConductance);
            state.Matrix[B1gNode, B1gNode] += (gcggb);
            state.Matrix[B1sNode, B1sNode] += (B1sourceConductance);
            state.Matrix[B1bNode, B1bNode] += (gbd + gbs - gcbgb - gcbdb - gcbsb);
            state.Matrix[B1dNodePrime, B1dNodePrime] +=
                (B1drainConductance + gds + gbd + xrev * (gm + gmbs) + gcddb);
            state.Matrix[B1sNodePrime, B1sNodePrime] +=
                (B1sourceConductance + gds + gbs + xnrm * (gm + gmbs) + gcssb);
            state.Matrix[B1dNode, B1dNodePrime] += (-B1drainConductance);
            state.Matrix[B1gNode, B1bNode] += (-gcggb - gcgdb - gcgsb);
            state.Matrix[B1gNode, B1dNodePrime] += (gcgdb);
            state.Matrix[B1gNode, B1sNodePrime] += (gcgsb);
            state.Matrix[B1sNode, B1sNodePrime] += (-B1sourceConductance);
            state.Matrix[B1bNode, B1gNode] += (gcbgb);
            state.Matrix[B1bNode, B1dNodePrime] += (-gbd + gcbdb);
            state.Matrix[B1bNode, B1sNodePrime] += (-gbs + gcbsb);
            state.Matrix[B1dNodePrime, B1dNode] += (-B1drainConductance);
            state.Matrix[B1dNodePrime, B1gNode] += ((xnrm - xrev) * gm + gcdgb);
            state.Matrix[B1dNodePrime, B1bNode] += (-gbd + (xnrm - xrev) * gmbs - gcdgb - gcddb - gcdsb);
            state.Matrix[B1dNodePrime, B1sNodePrime] += (-gds - xnrm * (gm + gmbs) + gcdsb);
            state.Matrix[B1sNodePrime, B1gNode] += (-(xnrm - xrev) * gm + gcsgb);
            state.Matrix[B1sNodePrime, B1sNode] += (-B1sourceConductance);
            state.Matrix[B1sNodePrime, B1bNode] += (-gbs - (xnrm - xrev) * gmbs - gcsgb - gcsdb - gcssb);
            state.Matrix[B1sNodePrime, B1dNodePrime] += (-gds - xrev * (gm + gmbs) + gcsdb);

        line1000: ;
        }

        /// <summary>
        /// Evaluate the B1 current
        /// </summary>
        /// <returns></returns>
        private void B1evaluate(Circuit ckt, BSIM1Model model, double vds, double vbs, double vgs, ref double gm, ref double gds, ref double gmbs,
            ref double qg, ref double qb, ref double qd, ref double cggb, ref double cgdb, ref double cgsb,
            ref double cbgb, ref double cbdb, ref double cbsb, ref double cdgb, ref double cddb,
            ref double cdsb, ref double cdrain, ref double von, ref double vdsat)
        {
            double Vfb;
            double Phi;
            double K1;
            double K2;
            double Vdd;
            double Ugs;
            double Uds;
            double dUgsdVbs;
            double Leff;
            double dUdsdVbs;
            double dUdsdVds;
            double Eta;
            double dEtadVds;
            double dEtadVbs;
            double Vpb;
            double SqrtVpb;
            double Von;
            double Vth;
            double dVthdVbs;
            double dVthdVds;
            double Vgs_Vth;
            double DrainCurrent;
            double G;
            double A;
            double Arg;
            double dGdVbs;
            double dAdVbs;
            double Beta;
            double Beta_Vds_0;
            double BetaVdd;
            double dBetaVdd_dVds;
            double Beta0;
            double dBeta0dVds;
            double dBeta0dVbs;
            double VddSquare;
            double C1;
            double C2;
            double dBetaVdd_dVbs;
            double dBeta_Vds_0_dVbs;
            double dC1dVbs;
            double dC2dVbs;
            double dBetadVgs;
            double dBetadVds;
            double dBetadVbs;
            double VdsSat = 0.0;
            double Argl1;
            double Argl2;
            double Vc;
            double Term1;
            double K;
            double Args1;
            double dVcdVgs;
            double dVcdVds;
            double dVcdVbs;
            double dKdVc;
            double dKdVgs;
            double dKdVds;
            double dKdVbs;
            double Args2;
            double Args3;
            double Warg1;
            double Vcut;
            double N;
            double N0;
            double NB;
            double ND;
            double Warg2;
            double Wds;
            double Wgs;
            double Ilimit;
            double Iexp;
            double Temp1;
            double Vth0;
            double Arg1;
            double Arg2;
            double Arg3;
            double Arg5;
            double Ent;
            double Vcom;
            double Vgb;
            double Vgb_Vfb;
            double VdsPinchoff;
            double EntSquare;
            double Vgs_VthSquare;
            double Argl3;
            double Argl4;
            double Argl5;
            double Argl6;
            double Argl7;
            double Argl8;
            double Argl9;
            double dEntdVds;
            double dEntdVbs;
            double cgbb;
            double cdbb;
            double cbbb;
            double WLCox;
            double Vtsquare;
            double Temp3;
            bool ChargeComputationNeeded = false;
            double co4v15;

            if ((ckt.Mode & (Circuit.Modes.Ac | Circuit.Modes.Tran)) != 0)
                ChargeComputationNeeded = true;
            else if ((ckt.Mode & Circuit.Modes.InitSmSig) != 0)
                ChargeComputationNeeded = true;
            else if (((ckt.Mode & Circuit.Modes.TranOp) != 0) && ((ckt.Mode & Circuit.Modes.Uic) != 0))
                ChargeComputationNeeded = true;
            else
                ChargeComputationNeeded = false;

            Vfb = B1vfb;
            Phi = B1phi;
            K1 = B1K1;
            K2 = B1K2;
            Vdd = model.B1vdd;
            if ((Ugs = B1ugs + B1ugsB * vbs) <= 0)
            {
                Ugs = 0;
                dUgsdVbs = 0.0;
            }
            else
            {
                dUgsdVbs = B1ugsB;
            }
            if ((Uds = B1uds + B1udsB * vbs +
                    B1udsD * (vds - Vdd)) <= 0)
            {
                Uds = 0.0;
                dUdsdVbs = dUdsdVds = 0.0;
            }
            else
            {
                Leff = B1l * 1.0e6 - model.B1deltaL; /* Leff in um */
                Uds = Uds / Leff;
                dUdsdVbs = B1udsB / Leff;
                dUdsdVds = B1udsD / Leff;
            }
            Eta = B1eta + B1etaB * vbs + B1etaD *
                (vds - Vdd);
            if (Eta <= 0)
            {
                Eta = 0;
                dEtadVds = dEtadVbs = 0.0;
            }
            else if (Eta > 1)
            {
                Eta = 1;
                dEtadVds = dEtadVbs = 0;
            }
            else
            {
                dEtadVds = B1etaD;
                dEtadVbs = B1etaB;
            }
            if (vbs < 0)
            {
                Vpb = Phi - vbs;
            }
            else
            {
                Vpb = Phi;
            }
            SqrtVpb = Math.Sqrt(Vpb);
            Von = Vfb + Phi + K1 * SqrtVpb - K2 * Vpb - Eta * vds;
            Vth = Von;
            dVthdVds = -Eta - dEtadVds * vds;
            dVthdVbs = K2 - 0.5 * K1 / SqrtVpb - dEtadVbs * vds;
            Vgs_Vth = vgs - Vth;

            G = 1.0 - 1.0 / (1.744 + 0.8364 * Vpb);
            A = 1.0 + 0.5 * G * K1 / SqrtVpb;
            A = Math.Max(A, 1.0);   /* Modified */
            Arg = Math.Max((1 + Ugs * Vgs_Vth), 1.0);
            dGdVbs = -0.8364 * (1 - G) * (1 - G);
            dAdVbs = 0.25 * K1 / SqrtVpb * (2 * dGdVbs + G / Vpb);

            if (Vgs_Vth < 0)
            {
                /* cutoff */
                DrainCurrent = 0;
                gm = 0;
                gds = 0;
                gmbs = 0;
                goto SubthresholdComputation;
            }

            /* Quadratic Interpolation for Beta0 (Beta at vgs  =  0, vds=Vds) */

            Beta_Vds_0 = (B1betaZero + B1betaZeroB * vbs);
            BetaVdd = (B1betaVdd + B1betaVddB * vbs);
            dBetaVdd_dVds = Math.Max(B1betaVddD, 0.0); /* Modified */
            if (vds > Vdd)
            {
                Beta0 = BetaVdd + dBetaVdd_dVds * (vds - Vdd);
                dBeta0dVds = dBetaVdd_dVds;
                dBeta0dVbs = B1betaVddB;
            }
            else
            {
                VddSquare = Vdd * Vdd;
                C1 = (-BetaVdd + Beta_Vds_0 + dBetaVdd_dVds * Vdd) / VddSquare;
                C2 = 2 * (BetaVdd - Beta_Vds_0) / Vdd - dBetaVdd_dVds;
                dBeta_Vds_0_dVbs = B1betaZeroB;
                dBetaVdd_dVbs = B1betaVddB;
                dC1dVbs = (dBeta_Vds_0_dVbs - dBetaVdd_dVbs) / VddSquare;
                dC2dVbs = dC1dVbs * (-2) * Vdd;
                Beta0 = (C1 * vds + C2) * vds + Beta_Vds_0;
                dBeta0dVds = 2 * C1 * vds + C2;
                dBeta0dVbs = dC1dVbs * vds * vds + dC2dVbs * vds + dBeta_Vds_0_dVbs;
            }

            /*Beta  =  Beta0 / ( 1 + Ugs * Vgs_Vth );*/

            Beta = Beta0 / Arg;
            dBetadVgs = -Beta * Ugs / Arg;
            dBetadVds = dBeta0dVds / Arg - dBetadVgs * dVthdVds;
            dBetadVbs = dBeta0dVbs / Arg + Beta * Ugs * dVthdVbs / Arg -
                Beta * Vgs_Vth * dUgsdVbs / Arg;

            /*VdsSat  = Math.Max( Vgs_Vth / ( A + Uds * Vgs_Vth ),  0.0);*/

            if ((Vc = Uds * Vgs_Vth / A) < 0.0) Vc = 0.0;
            Term1 = Math.Sqrt(1 + 2 * Vc);
            K = 0.5 * (1 + Vc + Term1);
            VdsSat = Math.Max(Vgs_Vth / (A * Math.Sqrt(K)), 0.0);

            if (vds < VdsSat)
            {
                /* Triode Region */
                /*Argl1  =  1 + Uds * vds;*/
                Argl1 = Math.Max((1 + Uds * vds), 1.0);
                Argl2 = Vgs_Vth - 0.5 * A * vds;
                DrainCurrent = Beta * Argl2 * vds / Argl1;
                gm = (dBetadVgs * Argl2 * vds + Beta * vds) / Argl1;
                gds = (dBetadVds * Argl2 * vds + Beta *
                    (Vgs_Vth - vds * dVthdVds - A * vds) -
                    DrainCurrent * (vds * dUdsdVds + Uds)) / Argl1;
                gmbs = (dBetadVbs * Argl2 * vds + Beta * vds *
                    (-dVthdVbs - 0.5 * vds * dAdVbs) -
                    DrainCurrent * vds * dUdsdVbs) / Argl1;
            }
            else
            {
                /* Pinchoff (Saturation) Region */
                Args1 = 1.0 + 1.0 / Term1;
                dVcdVgs = Uds / A;
                dVcdVds = Vgs_Vth * dUdsdVds / A - dVcdVgs * dVthdVds;
                dVcdVbs = (Vgs_Vth * dUdsdVbs - Uds *
                    (dVthdVbs + Vgs_Vth * dAdVbs / A)) / A;
                dKdVc = 0.5 * Args1;
                dKdVgs = dKdVc * dVcdVgs;
                dKdVds = dKdVc * dVcdVds;
                dKdVbs = dKdVc * dVcdVbs;
                Args2 = Vgs_Vth / A / K;
                Args3 = Args2 * Vgs_Vth;
                DrainCurrent = 0.5 * Beta * Args3;
                gm = 0.5 * Args3 * dBetadVgs + Beta * Args2 -
                    DrainCurrent * dKdVgs / K;
                gds = 0.5 * Args3 * dBetadVds - Beta * Args2 * dVthdVds -
                    DrainCurrent * dKdVds / K;
                gmbs = 0.5 * dBetadVbs * Args3 - Beta * Args2 * dVthdVbs -
                    DrainCurrent * (dAdVbs / A + dKdVbs / K);
            }

        SubthresholdComputation:

            N0 = B1subthSlope;
            Vcut = -40.0 * N0 * Circuit.CONSTvt0;

            /* The following 'if' statement has been modified so that subthreshold  *
             * current computation is always executed unless N0 >= 200. This should *
             * get rid of the Ids kink seen on Ids-Vgs plots at low Vds.            *
             *                                                Peter M. Lee          *
             *                                                6/8/90                *
             *  Old 'if' statement:                                                 *
             *  if( (N0 >=  200) || (Vgs_Vth < Vcut ) || (Vgs_Vth > (-0.5*Vcut)))   */

            if (N0 >= 200)
            {
                goto ChargeComputation;
            }

            NB = B1subthSlopeB;
            ND = B1subthSlopeD;
            N = N0 + NB * vbs + ND * vds; /* subthreshold slope */
            if (N < 0.5) N = 0.5;
            Warg1 = Math.Exp(-vds / Circuit.CONSTvt0);
            Wds = 1 - Warg1;
            Wgs = Math.Exp(Vgs_Vth / (N * Circuit.CONSTvt0));
            Vtsquare = Circuit.CONSTvt0 * Circuit.CONSTvt0;
            Warg2 = 6.04965 * Vtsquare * B1betaZero;
            Ilimit = 4.5 * Vtsquare * B1betaZero;
            Iexp = Warg2 * Wgs * Wds;
            DrainCurrent = DrainCurrent + Ilimit * Iexp / (Ilimit + Iexp);
            Temp1 = Ilimit / (Ilimit + Iexp);
            Temp1 = Temp1 * Temp1;
            Temp3 = Ilimit / (Ilimit + Wgs * Warg2);
            Temp3 = Temp3 * Temp3 * Warg2 * Wgs;
            /*    if ( Temp3 > Ilimit ) Temp3=Ilimit;*/
            gm = gm + Temp1 * Iexp / (N * Circuit.CONSTvt0);
            /* gds term has been modified to prevent blow up at Vds=0 */
            gds = gds + Temp3 * (-Wds / N / Circuit.CONSTvt0 * (dVthdVds +
                Vgs_Vth * ND / N) + Warg1 / Circuit.CONSTvt0);
            gmbs = gmbs - Temp1 * Iexp * (dVthdVbs + Vgs_Vth * NB / N) /
                (N * Circuit.CONSTvt0);

        ChargeComputation:

            /* Some Limiting of DC Parameters */
            if (DrainCurrent < 0.0) DrainCurrent = 0.0;
            if (gm < 0.0) gm = 0.0;
            if (gds < 0.0) gds = 0.0;
            if (gmbs < 0.0) gmbs = 0.0;

            WLCox = model.B1Cox *
                (B1l - model.B1deltaL * 1.0e-6) *
                (B1w - model.B1deltaW * 1.0e-6) * 1.0e4;   /* F */

            if (!ChargeComputationNeeded)
            {
                qg = 0;
                qd = 0;
                qb = 0;
                cggb = 0;
                cgsb = 0;
                cgdb = 0;
                cdgb = 0;
                cdsb = 0;
                cddb = 0;
                cbgb = 0;
                cbsb = 0;
                cbdb = 0;
                goto finished;
            }
            G = 1.0 - 1.0 / (1.744 + 0.8364 * Vpb);
            A = 1.0 + 0.5 * G * K1 / SqrtVpb;
            A = Math.Max(A, 1.0);   /* Modified */
            /*Arg  =  1 + Ugs * Vgs_Vth;*/
            dGdVbs = -0.8364 * (1 - G) * (1 - G);
            dAdVbs = 0.25 * K1 / SqrtVpb * (2 * dGdVbs + G / Vpb);
            Phi = Math.Max(0.1, Phi);

            if (model.B1channelChargePartitionFlag >= 1)
            {

                /*0/100 partitioning for drain/source chArges at the saturation region*/
                Vth0 = Vfb + Phi + K1 * SqrtVpb;
                Vgs_Vth = vgs - Vth0;
                Arg1 = A * vds;
                Arg2 = Vgs_Vth - 0.5 * Arg1;
                Arg3 = vds - Arg1;
                Arg5 = Arg1 * Arg1;
                dVthdVbs = -0.5 * K1 / SqrtVpb;
                dAdVbs = 0.5 * K1 * (0.5 * G / Vpb - 0.8364 * (1 - G) * (1 - G)) /
                    SqrtVpb;
                Ent = Math.Max(Arg2, 1.0e-8);
                dEntdVds = -0.5 * A;
                dEntdVbs = -dVthdVbs - 0.5 * vds * dAdVbs;
                Vcom = Vgs_Vth * Vgs_Vth / 6.0 - 1.25e-1 * Arg1 *
                    Vgs_Vth + 2.5e-2 * Arg5;
                VdsPinchoff = Math.Max(Vgs_Vth / A, 0.0);
                Vgb = vgs - vbs;
                Vgb_Vfb = Vgb - Vfb;

                if (Vgb_Vfb < 0)
                {
                    /* Accumulation Region */
                    qg = WLCox * Vgb_Vfb;
                    qb = -qg;
                    qd = 0.0;
                    cggb = WLCox;
                    cgdb = 0.0;
                    cgsb = 0.0;
                    cbgb = -WLCox;
                    cbdb = 0.0;
                    cbsb = 0.0;
                    cdgb = 0.0;
                    cddb = 0.0;
                    cdsb = 0.0;
                    goto finished;
                }
                else if (vgs < Vth0)
                {
                    /* Subthreshold Region */
                    qg = 0.5 * WLCox * K1 * K1 * (-1 +
                        Math.Sqrt(1 + 4 * Vgb_Vfb / (K1 * K1)));
                    qb = -qg;
                    qd = 0.0;
                    cggb = WLCox / Math.Sqrt(1 + 4 * Vgb_Vfb / (K1 * K1));
                    cgdb = cgsb = 0.0;
                    cbgb = -cggb;
                    cbdb = cbsb = cdgb = cddb = cdsb = 0.0;
                    goto finished;
                }
                else if (vds < VdsPinchoff)
                {    /* triode region  */
                    /*Vgs_Vth2 = Vgs_Vth*Vgs_Vth;*/
                    EntSquare = Ent * Ent;
                    Argl1 = 1.2e1 * EntSquare;
                    Argl2 = 1.0 - A;
                    Argl3 = Arg1 * vds;
                    /*Argl4 = Vcom/Ent/EntSquare;*/
                    if (Ent > 1.0e-8)
                    {
                        Argl5 = Arg1 / Ent;
                        /*Argl6 = Vcom/EntSquare;*/
                    }
                    else
                    {
                        Argl5 = 2.0;
                        Argl6 = 4.0 / 1.5e1;
                    }
                    Argl7 = Argl5 / 1.2e1;
                    Argl8 = 6.0 * Ent;
                    Argl9 = 0.125 * Argl5 * Argl5;
                    qg = WLCox * (vgs - Vfb - Phi - 0.5 * vds + vds * Argl7);
                    qb = WLCox * (-Vth0 + Vfb + Phi + 0.5 * Arg3 - Arg3 * Argl7);
                    qd = -WLCox * (0.5 * Vgs_Vth - 0.75 * Arg1 +
                        0.125 * Arg1 * Argl5);
                    cggb = WLCox * (1.0 - Argl3 / Argl1);
                    cgdb = WLCox * (-0.5 + Arg1 / Argl8 - Argl3 * dEntdVds /
                        Argl1);
                    cgbb = WLCox * (vds * vds * dAdVbs * Ent - Argl3 * dEntdVbs) /
                        Argl1;
                    cgsb = -(cggb + cgdb + cgbb);
                    cbgb = WLCox * Argl3 * Argl2 / Argl1;
                    cbdb = WLCox * Argl2 * (0.5 - Arg1 / Argl8 + Argl3 * dEntdVds /
                        Argl1);
                    cbbb = -WLCox * (dVthdVbs + 0.5 * vds * dAdVbs + vds *
                        vds * ((1.0 - 2.0 * A) * dAdVbs * Ent - Argl2 *
                        A * dEntdVbs) / Argl1);
                    cbsb = -(cbgb + cbdb + cbbb);
                    cdgb = -WLCox * (0.5 - Argl9);
                    cddb = WLCox * (0.75 * A - 0.25 * A * Arg1 / Ent +
                        Argl9 * dEntdVds);
                    cdbb = WLCox * (0.5 * dVthdVbs + vds * dAdVbs *
                        (0.75 - 0.25 * Argl5) + Argl9 * dEntdVbs);
                    cdsb = -(cdgb + cddb + cdbb);
                    goto finished;
                }
                else if (vds >= VdsPinchoff)
                {    /* saturation region   */
                    Args1 = 1.0 / (3.0 * A);
                    qg = WLCox * (vgs - Vfb - Phi - Vgs_Vth * Args1);
                    qb = WLCox * (Vfb + Phi - Vth0 + (1.0 - A) * Vgs_Vth * Args1);
                    qd = 0.0;
                    cggb = WLCox * (1.0 - Args1);
                    cgdb = 0.0;
                    cgbb = WLCox * Args1 * (dVthdVbs + Vgs_Vth * dAdVbs / A);
                    cgsb = -(cggb + cgdb + cgbb);
                    cbgb = WLCox * (Args1 - 1.0 / 3.0);
                    cbdb = 0.0;
                    cbbb = -WLCox * ((2.0 / 3.0 + Args1) * dVthdVbs +
                        Vgs_Vth * Args1 * dAdVbs / A);      /* Modified */
                    cbsb = -(cbgb + cbdb + cbbb);
                    cdgb = 0.0;
                    cddb = 0.0;
                    cdsb = 0.0;
                    goto finished;
                }

                goto finished;

            }
            else
            {
                /* ChannelChargePartionFlag  < = 0 */

                /*40/60 partitioning for drain/source chArges at the saturation region*/
                co4v15 = 4.0 / 15.0;
                Vth0 = Vfb + Phi + K1 * SqrtVpb;
                Vgs_Vth = vgs - Vth0;
                Arg1 = A * vds;
                Arg2 = Vgs_Vth - 0.5 * Arg1;
                Arg3 = vds - Arg1;
                Arg5 = Arg1 * Arg1;
                dVthdVbs = -0.5 * K1 / SqrtVpb;
                dAdVbs = 0.5 * K1 * (0.5 * G / Vpb - 0.8364 * (1 - G) * (1 - G)) / SqrtVpb;
                Ent = Math.Max(Arg2, 1.0e-8);
                dEntdVds = -0.5 * A;
                dEntdVbs = -dVthdVbs - 0.5 * vds * dAdVbs;
                Vcom = Vgs_Vth * Vgs_Vth / 6.0 - 1.25e-1 * Arg1 * Vgs_Vth + 2.5e-2 * Arg5;
                VdsPinchoff = Math.Max(Vgs_Vth / A, 0.0);
                Vgb = vgs - vbs;
                Vgb_Vfb = Vgb - Vfb;

                if (Vgb_Vfb < 0)
                {           /* Accumulation Region */
                    qg = WLCox * Vgb_Vfb;
                    qb = -qg;
                    qd = 0.0;
                    cggb = WLCox;
                    cgdb = 0.0;
                    cgsb = 0.0;
                    cbgb = -WLCox;
                    cbdb = 0.0;
                    cbsb = 0.0;
                    cdgb = 0.0;
                    cddb = 0.0;
                    cdsb = 0.0;
                    goto finished;
                }
                else if (vgs < Vth0)
                {    /* Subthreshold Region */
                    qg = 0.5 * WLCox * K1 * K1 * (-1 + Math.Sqrt(1 + 4 * Vgb_Vfb / (K1 * K1)));
                    qb = -qg;
                    qd = 0.0;
                    cggb = WLCox / Math.Sqrt(1 + 4 * Vgb_Vfb / (K1 * K1));
                    cgdb = cgsb = 0.0;
                    cbgb = -cggb;
                    cbdb = cbsb = cdgb = cddb = cdsb = 0.0;
                    goto finished;
                }
                else if (vds < VdsPinchoff)
                {      /* triode region */

                    Vgs_VthSquare = Vgs_Vth * Vgs_Vth;
                    EntSquare = Ent * Ent;
                    Argl1 = 1.2e1 * EntSquare;
                    Argl2 = 1.0 - A;
                    Argl3 = Arg1 * vds;
                    Argl4 = Vcom / Ent / EntSquare;
                    if (Ent > 1.0e-8)
                    {
                        Argl5 = Arg1 / Ent;
                        Argl6 = Vcom / EntSquare;
                    }
                    else
                    {
                        Argl5 = 2.0;
                        Argl6 = 4.0 / 1.5e1;
                    }
                    Argl7 = Argl5 / 1.2e1;
                    Argl8 = 6.0 * Ent;
                    qg = WLCox * (vgs - Vfb - Phi - 0.5 * vds + vds * Argl7);
                    qb = WLCox * (-Vth0 + Vfb + Phi + 0.5 * Arg3 - Arg3 * Argl7);
                    qd = -WLCox * (0.5 * (Vgs_Vth - Arg1) + Arg1 * Argl6);
                    cggb = WLCox * (1.0 - Argl3 / Argl1);
                    cgdb = WLCox * (-0.5 + Arg1 / Argl8 - Argl3 * dEntdVds / Argl1);
                    cgbb = WLCox * (vds * vds * dAdVbs * Ent - Argl3 * dEntdVbs) / Argl1;
                    cgsb = -(cggb + cgdb + cgbb);
                    cbgb = WLCox * Argl3 * Argl2 / Argl1;
                    cbdb = WLCox * Argl2 * (0.5 - Arg1 / Argl8 + Argl3 * dEntdVds / Argl1);
                    cbbb = -WLCox * (dVthdVbs + 0.5 * vds * dAdVbs + vds * vds * ((1.0 - 2.0 * A)
                        * dAdVbs * Ent - Argl2 * A * dEntdVbs) / Argl1);
                    cbsb = -(cbgb + cbdb + cbbb);
                    cdgb = -WLCox * (0.5 + Arg1 * (4.0 * Vgs_Vth - 1.5 * Arg1) / Argl1 -
                        2.0 * Arg1 * Argl4);
                    cddb = WLCox * (0.5 * A + 2.0 * Arg1 * dEntdVds * Argl4 - A * (2.0 * Vgs_VthSquare
                        - 3.0 * Arg1 * Vgs_Vth + 0.9 * Arg5) / Argl1);
                    cdbb = WLCox * (0.5 * dVthdVbs + 0.5 * vds * dAdVbs + 2.0 * Arg1 * dEntdVbs
                        * Argl4 - vds * (2.0 * Vgs_VthSquare * dAdVbs - 4.0 * A * Vgs_Vth * dVthdVbs - 3.0
                        * Arg1 * Vgs_Vth * dAdVbs + 1.5 * A * Arg1 * dVthdVbs + 0.9 * Arg5 * dAdVbs)
                        / Argl1);
                    cdsb = -(cdgb + cddb + cdbb);
                    goto finished;
                }
                else if (vds >= VdsPinchoff)
                {      /* saturation region */

                    Args1 = 1.0 / (3.0 * A);
                    qg = WLCox * (vgs - Vfb - Phi - Vgs_Vth * Args1);
                    qb = WLCox * (Vfb + Phi - Vth0 + (1.0 - A) * Vgs_Vth * Args1);
                    qd = -co4v15 * WLCox * Vgs_Vth;
                    cggb = WLCox * (1.0 - Args1);
                    cgdb = 0.0;
                    cgbb = WLCox * Args1 * (dVthdVbs + Vgs_Vth * dAdVbs / A);
                    cgsb = -(cggb + cgdb + cgbb);
                    cbgb = WLCox * (Args1 - 1.0 / 3.0);
                    cbdb = 0.0;
                    cbbb = -WLCox * ((2.0 / 3.0 + Args1) * dVthdVbs + Vgs_Vth * Args1 * dAdVbs / A);
                    cbsb = -(cbgb + cbdb + cbbb);
                    cdgb = -co4v15 * WLCox;
                    cddb = 0.0;
                    cdbb = co4v15 * WLCox * dVthdVbs;
                    cdsb = -(cdgb + cddb + cdbb);
                    goto finished;
                }
            }

        finished:       /* returning Values to Calling Routine */
            gm = Math.Max(gm, 0.0);
            gds = Math.Max(gds, 0.0);
            gmbs = Math.Max(gmbs, 0.0);
            cdrain = Math.Max(DrainCurrent, 0.0);
            von = Von;
            vdsat = VdsSat;
            // vdsat = VdsSat; // NOTE: "Use of unassigned local variable VdsSat"
        }

        /// <summary>
        /// Evaluate the B1 capacitors
        /// </summary>
        private void B1mosCap(Circuit ckt, double vgd, double vgs, double vgb, double[] args,
            double cbgb, double cbdb, double cbsb, double cdgb, double cddb, double cdsb,
            ref double gcggb, ref double gcgdb, ref double gcgsb, ref double gcbgb, ref double gcbdb,
            ref double gcbsb, ref double gcdgb, ref double gcddb, ref double gcdsb,
            ref double gcsgb, ref double gcsdb, ref double gcssb, ref double qGate, ref double qBulk, ref double qDrain, ref double qSource)
        {
            double ag0 = ckt.Integration.Deriv; // ckt.CKTag[0];

            /* compute equivalent conductance */
            gcdgb = (cdgb - args[0]) * ag0;
            gcddb = (cddb + args[3] + args[0]) * ag0;
            gcdsb = cdsb * ag0;
            gcsgb = -(args[5] + cbgb + cdgb + args[1]) * ag0;
            gcsdb = -(args[6] + cbdb + cddb) * ag0;
            gcssb = (args[4] + args[1] -
                (args[7] + cbsb + cdsb)) * ag0;
            gcggb = (args[5] + args[0] +
                args[1] + args[2]) * ag0;
            gcgdb = (args[6] - args[0]) * ag0;
            gcgsb = (args[7] - args[1]) * ag0;
            gcbgb = (cbgb - args[2]) * ag0;
            gcbdb = (cbdb - args[3]) * ag0;
            gcbsb = (cbsb - args[4]) * ag0;

            /* compute total terminal charge */
            double qgd = args[0] * vgd;
            double qgs = args[1] * vgs;
            double qgb = args[2] * vgb;
            qGate = qGate + qgd + qgs + qgb;
            qBulk = qBulk - qgb;
            qDrain = qDrain - qgd;
            qSource = -(qGate + qBulk + qDrain);
        }

        /// <summary>
        /// Load the matrix with the BSIM1 model for AC or Pole-Zero analysis
        /// </summary>
        /// <param name="ckt"></param>
        /// <param name="p"></param>
        public override void AcLoad(Circuit ckt, Model model, object p = null)
        {
            CircuitStateComplex state = ckt.StateComplex;
            CircuitState cstate = ckt.State;

            int xnrm;
            int xrev;
            double gdpr;
            double gspr;
            double gm;
            double gds;
            double gmbs;
            double gbd;
            double gbs;
            double capbd;
            double capbs;
            double xcggb;
            double xcgdb;
            double xcgsb;
            double xcbgb;
            double xcbdb;
            double xcbsb;
            double xcddb;
            double xcssb;
            double xcdgb;
            double xcsgb;
            double xcdsb;
            double xcsdb;
            double cggb;
            double cgdb;
            double cgsb;
            double cbgb;
            double cbdb;
            double cbsb;
            double cddb;
            double cdgb;
            double cdsb;

            if (B1mode >= 0)
            {
                xnrm = 1;
                xrev = 0;
            }
            else
            {
                xnrm = 0;
                xrev = 1;
            }
            gdpr = B1drainConductance;
            gspr = B1sourceConductance;
            gm = cstate.States[0][B1state + B1gm];
            gds = cstate.States[0][B1state + B1gds];
            gmbs = cstate.States[0][B1state + B1gmbs];
            gbd = cstate.States[0][B1state + B1gbd];
            gbs = cstate.States[0][B1state + B1gbs];
            capbd = cstate.States[0][B1state + B1capbd];
            capbs = cstate.States[0][B1state + B1capbs];
            /*
             *    charge oriented model parameters
             */

            cggb = cstate.States[0][B1state + B1cggb];
            cgsb = cstate.States[0][B1state + B1cgsb];
            cgdb = cstate.States[0][B1state + B1cgdb];

            cbgb = cstate.States[0][B1state + B1cbgb];
            cbsb = cstate.States[0][B1state + B1cbsb];
            cbdb = cstate.States[0][B1state + B1cbdb];

            cdgb = cstate.States[0][B1state + B1cdgb];
            cdsb = cstate.States[0][B1state + B1cdsb];
            cddb = cstate.States[0][B1state + B1cddb];

            xcdgb = (cdgb - B1GDoverlapCap);
            xcddb = (cddb + capbd + B1GDoverlapCap);
            xcdsb = cdsb;
            xcsgb = -(cggb + cbgb + cdgb + B1GSoverlapCap);
            xcsdb = -(cgdb + cbdb + cddb);
            xcssb = (capbs + B1GSoverlapCap - (cgsb + cbsb + cdsb));
            xcggb = (cggb + B1GDoverlapCap + B1GSoverlapCap +
                    B1GBoverlapCap);
            xcgdb = (cgdb - B1GDoverlapCap);
            xcgsb = (cgsb - B1GSoverlapCap);
            xcbgb = (cbgb - B1GBoverlapCap);
            xcbdb = (cbdb - capbd);
            xcbsb = (cbsb - capbs);

            state.Matrix[B1gNode, B1gNode] += xcggb * state.Omega;
            state.Matrix[B1bNode, B1bNode] += (-xcbgb - xcbdb - xcbsb) * state.Omega;
            state.Matrix[B1dNodePrime, B1dNodePrime] += xcddb * state.Omega;
            state.Matrix[B1sNodePrime, B1sNodePrime] += xcssb * state.Omega;
            state.Matrix[B1gNode, B1bNode] += (-xcggb - xcgdb - xcgsb) * state.Omega;
            state.Matrix[B1gNode, B1dNodePrime] += xcgdb * state.Omega;
            state.Matrix[B1gNode, B1sNodePrime] += xcgsb * state.Omega;
            state.Matrix[B1bNode, B1gNode] += xcbgb * state.Omega;
            state.Matrix[B1bNode, B1dNodePrime] += xcbdb * state.Omega;
            state.Matrix[B1bNode, B1sNodePrime] += xcbsb * state.Omega;
            state.Matrix[B1dNodePrime, B1gNode] += xcdgb * state.Omega;
            state.Matrix[B1dNodePrime, B1bNode] += (-xcdgb - xcddb - xcdsb) * state.Omega;
            state.Matrix[B1dNodePrime, B1sNodePrime] += xcdsb * state.Omega;
            state.Matrix[B1sNodePrime, B1gNode] += xcsgb * state.Omega;
            state.Matrix[B1sNodePrime, B1bNode] += (-xcsgb - xcsdb - xcssb) * state.Omega;
            state.Matrix[B1sNodePrime, B1dNodePrime] += xcsdb * state.Omega;
            state.Matrix[B1dNode, B1dNode] += gdpr;
            state.Matrix[B1sNode, B1sNode] += gspr;
            state.Matrix[B1bNode, B1bNode] += gbd + gbs;
            state.Matrix[B1dNodePrime, B1dNodePrime] += gdpr + gds + gbd + xrev * (gm + gmbs);
            state.Matrix[B1sNodePrime, B1sNodePrime] += gspr + gds + gbs + xnrm * (gm + gmbs);
            state.Matrix[B1dNode, B1dNodePrime] -= gdpr;
            state.Matrix[B1sNode, B1sNodePrime] -= gspr;
            state.Matrix[B1bNode, B1dNodePrime] -= gbd;
            state.Matrix[B1bNode, B1sNodePrime] -= gbs;
            state.Matrix[B1dNodePrime, B1dNode] -= gdpr;
            state.Matrix[B1dNodePrime, B1gNode] += (xnrm - xrev) * gm;
            state.Matrix[B1dNodePrime, B1bNode] += -gbd + (xnrm - xrev) * gmbs;
            state.Matrix[B1dNodePrime, B1sNodePrime] += -gds - xnrm * (gm + gmbs);
            state.Matrix[B1sNodePrime, B1gNode] += -(xnrm - xrev) * gm;
            state.Matrix[B1sNodePrime, B1sNode] -= gspr;
            state.Matrix[B1sNodePrime, B1bNode] += -gbs - (xnrm - xrev) * gmbs;
            state.Matrix[B1sNodePrime, B1dNodePrime] += -gds - xrev * (gm + gmbs);
        }
    }
}
